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Abstract 

We discuss modern numerical methods for quantum spin systems and their 
apphcation to magnetic molecules. 
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1. Introduction 

The knowledge of the energy eigenvalues and eigenstates of small mag- 
netic systems such as magnetic molecules is indispensable for a complete 
understanding of their spectroscopic, dynamic, and thermodynamic proper- 
ties. In this respect the numerical exact diagonalization of the appropriate 
quantum Hamiltonian is the method of choice. Nevertheless, such an attempt 
is very often severely restricted due to the huge dimension of the underlying 
Hilbert space. For a magnetic system of N spins of spin quantum number s 
the dimension is (2s -|- 1)^ which grows exponentially with N . 

Group theoretical methods can help to case this numerical problem. 
Along these lines much effort has been put into the development of an efficient 
numerical diagonalization technique of the Heisenbcrg model 



using irreducible tensor operators and thus employing SU(2) symmetry of 
angular momenta [1, 2, 3, 4, 5, 6]. A combination of this meanwhile well es- 
tablished technique with point-group symmetries could be developed over the 
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past years, first for those point-group symmetries that are compatible with 
the spin couphng scheme, i.e. avoid comphcated basis transforms between 
different couphng schemes [7, 8, 9, 10], later also for general point groups 
[11, 12, 13]. Nevertheless, if the dimension of the largest Hilbert subspace 
exceeds about 10^ complete numerically exact diagonalization is no longer 
possible with current computers and programs. 

Fortunately, a few very accurate approximations have been developed 
that can be applied to quantum spin systems. In this article we are going 
to discuss the Finite Temperature Lanczos Method (FTLM), the Density 
Matrix Renormalization Group (DMRG) and its dynamical variant as well 
as Quantum Monte Carlo (QMC). 

For problems with Hilbert space dimensions of up to roughly lO^'' - the 
Finite Temperature Lanczos Method (FTLM) provides a very accurate and 
astonishingly easy to program method [14, 15]. In recent publications we 
demonstrated that this method is indeed capable of evaluating thermody- 
namic observables for magnetic molecules with an accuracy that is nearly 
indistinguishable from exact results [16, 17, 18]. So far we encountered only 
one problem where achieving a satisfying accuracy posed a problem - the 
Fe™ ferric wheel [19]. Mathematically this methods relies on the idea of 
trace estimators [20] ; it is not restricted to spin systems and for instance also 
applied in quantum chemistry [21, 22]. 

The Density Matrix Renormalization Group (DMRG) method is a vari- 
ational method that approximates the true eigenstates by so-called matrix- 
product states [23, 24, 25, 26]. These states are iteratively constructed, thus 
the method allows to treat the full Heisenberg Hamiltonian but in a reduced 
Hilbert space. The Hilbert space is truncated in a controlled way and the 
accuracy of the method can be estimated with the help of a truncation error. 
In the field of molecular magnetism DMRG has been applied for instance to 
the Heisenberg icosidodecahedron with s = 5/2 [27, 28], i.e. a model of the 
Mo72Fe3o Keplerate [29, 30]. 

DMRG can be extended in order to evaluate transition matrix elements. 
This variant is called Dynamical DMRG (DDMRG); it aims at the calculation 
of dynamical correlation functions as needed for the description of Inelastic 
Neutron Scattering cross sections [31, 32]. In a recent article we could show 
that DDMRG is able to model INS spectra of the very large magnetic ring 
molecule Fcig with unprecedented accuracy and thus allows to determine 
model parameters which would be impossible using only observables such as 
susceptibility [33]. This extension is not discussed in this article. 
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Finally we would like to provide examples for the application of Quantum 
Monte Carlo (QMC) [34, 35, 36] to magnetic molecules. This approximate 
method which works accurately only for non-frustrated quantum spin sys- 
tems [37], has already been applied to several molecular systems by Larry 
Engelhardt [38, 39, 40, 41]. He also provides a very popular program - FIT- 
MART - with which one can deduce Heisenberg exchange parameters from 
susceptibility data [42]. 

Some of the discussed methods are freely available as program packages. 
Besides FIT-MART the program MAGPACK [43] can be used for complete 
diagonalization. Approximate methods such as DMRG and QMC are pro- 
vided by the ALPS package [44, 45, 46]. 

The article is organized as follows. In Section 2 the Finite- Temperature 
Lanczos Method is introduced. Section 3 discusses the application of Quan- 
tum Monte Carlo, and Section 4 briefly introduces to the Density Matrix 
Renormalization Group. 

2. Application of the Finite-Temperature Lanczos Method to giant 
gadolinium clusters 




Figure 1: The core structure of {Gdi2Mo4} is a truncated tetrahedron. The bullets rep- 
resent the 12 spin sites and the edges correspond to the 18 exchange interactions between 
nearest-neighbor spins. The exchange inside the four triangles is named Ji, between tri- 
angles J2. 

In a recent publication heterometallic cluster complexes {Lni2Mo4} fea- 
turing a Lni2 core that has the structure of a distorted truncated tetrahe- 
dron, see Fig. 1, were reported [47]. The experimental magnetic studies of 
the {Gdi2Mo4} were accompanied on the theoretical side by calculations that 



3 



replaced the Gd spin of s = 7/2 by fictitious spins s = 5/2 since otherwise 
the calculation would not have been feasible in a reasonable time (of several 
weeks [sic!]). Now, after a few months, the calculations using the Finite- 
Temperature Lanczos Method for iV = 12 spins s = 7/2 are completed. 
Before presenting the results a short reminder of the method shall be given 
which in detail is explained elsewhere [14, 15, 16]. 

For the evaluation of thermodynamic properties in the canonical ensemble 
the exact partition function Z depending on temperature T and magnetic 
field B is given by 

= ^(^|e-^^^|^). (2) 

Here { denotes an orthonormal basis of the respective Hilbert space. 
Following the ideas of Refs. [14, 15] the unknown matrix elements are ap- 
proximated as 

Nl 
n=l 

which yields for the partition function 

z{T,E) » !^i:E'="''""i("Mi''>r w 

v=\ n=l 

For this procedure | ) is taken as the initial vector of a Lanczos iteration. 
This iteration consists of A^^ Lanczos steps, which span a respective Krylow 
space, in which the Hamiltonian is diagonalized. This yields the A^l Lanczos 
eigenvectors | n(v) ) as well as the associated Lanczos energy eigenvalues e^n^ . 
They are enumerated by n = 1, . . . , A''^. The number of Lanczos steps Ni, 
is a parameter of the approximation; ~ 100 is usually a good value. 
In addition, the complete and thus very large sum over all states \v) \^ 
replaced by a summation over a subset of R random vectors, where R is the 
second parameter of the method. For many cases R can be rather small, e.g. 
R ^ 20, whereas for other systems convergence is achieved only for R ^ 100. 
An observable would then be calculated as 

^ ' ^ r ^ 1^=1 n=l 

x{n{u,r)\o\u,r){u,r\n{u,r)) . (5) 
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Here F labels the irreducible representations of a symmetry group that can 
be used to split the Hilbert space into subspaces T-LiT) in order to increase 
the accuracy. In the following calculations we decomposed the Hilbert space 
according to the total magnetic quantum number M. 
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Figure 2: Magnetization of {Gdi2Mo4} as function of temperature (l.h.s.) as well as 
of applied magnetic field (r.h.s.). The experimental data [47] are given as symbols, the 
theoretical calculations as curves. 

The magnetization of {Gdi2Mo4} was evaluated for four different parame- 
ter sets. Since the total dimension is a staggering (2s+l)^ = 68, 719, 476, 736 
and even the dimension of the largest Hilbert subspace with M = is still 
3, 409, 213, 016, the calculations needed about a quarter of a year on a super- 
computer. As Fig. 2 shows, the exchange interactions are antiferromagnetic 
and of the order of -0.05 K. Since they are so small, the experimental data, 
taken from [47], is not sufficient to disentangle between scenarios where the 
interactions Ji between spins within triangles and J2 between triangles are 
the same or different. A scenario where only interactions between triangles 
bind the spins into dimers can be excluded, but a scenario where the system 
would consist of uncoupled triangles cannot be excluded. 

Figure 3 displays the magnetocaloric behavior for a coupling scheme with 
Ji = J2 = —0.05 K. The l.h.s. shows a set of isentropes, i.e. curves which 
the system would follow when the magnetic field is reduced in an adiabatic 
process. The figure on the r.h.s. shows the isothermal entropy changes for 
field sweeps from B = (1,2,7) T, respectively, to B = 0. The entropy 
differences are rather large at low temperature as expected for a weakly 
couple gadolinium system. 
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Figure 3: Theoretical niagnetic isentropes {Gdi2Mo4} as function of temperature and 
field (l.h.s.) and isothermal entropy change for various starting fields (r.h.s.). The chosen 
parameters are Ji = J 2 = —0.05 K. 



3. Quantum Monte Carlo 

Quantum Monte Carlo (QMC) [34, 35, 36] is a very powerful method 
for non-frustrated, i.e. bipartite quantum spin systems. For a discussion of 
frustration see e.g. [48]. The method can easily deal with up to 100 or more 
spins. In the field of molecule-based magnetism is was applied to several spin 
systems, e.g. homo- and heterometallic rings [38, 39, 40, 41] as well as to a 
one-dimensional spin tube [49]. In the latter publication the heat capacity 
of a system of = 100 spins with s = 3/2 was calculated with the help of 
QMC. 

Again the idea is to approximate the partition function. This time the 
partition function is chopped (sliced) in the sense that the exponential is 
written as a product of m exponentials with exponents divided by m (Trotter- 
Suzuki decomposition [50, 51, 52]). For m 00 the exponential can be 
written as a product of the exponentials of the (even non-commuting) parts 
of the Hamiltonian. One can as well linearize the exponential for large enough 
m. In any case, the multi-index sum is evaluated in a Monte-Carlo fashion 
as sketched in the equations below: 

ZiT,B) = 5^(z/|e-''5|z/)=5^(z/| [exp|-/3ff/m}]™|z/) (6) 
— 5^ ( I exp|— /3ff/m| |a)(Q;| exp|— /3i^/m| |7)(7| ••• 

v,a,7,... 
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As an example the magnetic susceptibility of the ten-membered ferric wheel, 
which was synthesized 18 years ago [53], is presented in Fig. 4. The symbols 
mark the experimental values [53], the solid curve shows the result of exact 
diagonalization [13], and the dotted curve shows the result obtained with 
ALPS QMC using 10^ steps for equilibration and lO^'' steps for the Monte- 
Carlo sampling for each temperature. As one can see the QMC result is 
indistinguishable from the exact one. 
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Figure 4: Susceptibility of an antiferromagnetically coupled spin ring with = 10 and 
s = 5/2. The exchange parameter J = —9.6 cm""'^ as well as the susceptibility data are 
taken from Ref. [53] . The solid curve displays the result of complete matrix diagonalization 
[13] whereas the dotted curve depicts the QMC result. 



4. DMRG results 

DMRG is a technically rather involved method. One should however keep 
in mind that it is a variational method that constructs trial states in a certain 
way. Being variational means that a lower energy for the trial ground state 
corresponds to having obtained a better approximation since the energy is 
bounded from below by the true value (Ritz's variational principle). A maybe 
helpful (mis-) conception of the idea of the method could be the following: 
Assume that you want to describe an eight-membered spin chain as depicted 
in Fig. 5. One starts with a small subsystem of two spins, diagonalizes the 
Hamiltonian and keeps only the lowest m eigenstates. Then one adds spins 
sequentially and every time sets up a new basis built of the old kept states 
and the states representing the added spin, diagonalizes the Hamiltonian and 
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Figure 5: Scheme of DMRG: first the system is iteratively enlarged up to its fuh size, then 
a sweep algorithm improves the accuracy further. 

again keeps the lowest m states. This idea is brilhant (and called Numerical 
Renormalization Group), except that it does not work in this naive fashion. 
Steve White found out, that instead of keeping the lowest m eigenstates 
of the Hamiltonian it is much better to keep m eigenstates of a reduced 
density matrix in order to iteratively built up the system and to represent 
the Hamiltonian [23]. The density matrix is given by the target state, e.g. 
the ground state |\l/),asp = |\E')(\I^| and to reduce it means to trace over a 
part of the system, i.e. one of the colored parts in Fig. 5. The representation 
can be further improved by running a so-called sweep algorithm in which 
the system is subdivided into unequal blocks for which the Hamiltonian is 
diagonalized and the density matrix calculated. 

For the non-expert this seems to be rather obscure, but contrary to several 
other methods, DMRG is (1) variational, (2) a controlled approximation, 
i.e. with m — 7- oo one approaches the exact result, and (3) offers accuracy 
estimators in the form of the truncated weight or the entanglement entropy. 
An extrapolation to the exact result is thus possible by using these measures. 
Although DMRG is best suited for open one-dimensional chain systems it 
can be applied to finite-size clusters, too. The resulting convergence, which 
is exponential in m for one-dimensional chains, is somewhat slower, e.g. like 
1/m for a spin cluster such as the icosidodecahedron [27]. 

Here we would like to demonstrate the power of the method by showing 
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theoretical magnetization curves for another ferric wheel, this time an Feig 
ring of = 18 spins with s = 5/2 [54]. The dimension of the Hilbert space for 
this system is about 10^^ which again renders an exact treatment impossible. 
The molecule was investigated by means of DMRG and DDMRG [33], the 
later results were compared to INS data and utilized to fix the parameters 
of the model. 
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Figure 6: Comparison of the zero-temperature magnetization curves for three different 
parameter sets as obtained by standard DMRG calculations for the N = 18, s ~ 5/2 spin 
ring. We have kept up to 3000 density matrix eigenstates for the Ji — J2 model and up 
to 1200 for the Ji ^ J2 models. The truncated weights are smaller than 10~®. 

In Fig. 6 we present the (T = 0) magnetization curves for the three 
parameterizations discussed in [33]. The results are very interesting. The 
magnetization curves are virtually identical up to 25 Tesla. The step widths 
of the magnetization curve for a Cis-symmetric ring (single J) are approx- 
imately the same for all steps (apart from the very last steps), as would 
also be the case for the rotational band approximation [55] . The magnetiza- 
tion curves for the other two models (Ce-symmetric: J1-J1-J2) deviate from 
this behavior. Up to approximately 100 T, the two models with Ji 7^ J2 
give very similar magnetization curves and considerable differences appear 
only at even higher fields. The magnetization curve for Ji = 1.42 meV and 
J2 = 4.57 meV shows two magnetization plateaus at higher fields. Plateaus 
in zero-temperature magnetization curves usually emerge in geometrically 
frustrated spin systems [56] . This system is, however, not geometrically frus- 
trated so that the emergence of a plateau is an interesting effect. 
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Figure 7: Dependence of the ^-component of the local magnetization (which is proportional 
to —{sj}) for two different parameter sets as a function of an external magnetic field 

B — Bcz- The results were obtained using standard DMRG and T = 0. The result 
for the uniform model is not shown because for this model the local moments are simply 
proportional to the magnetization curve presented in Fig. 6. "Symmetric position" denotes 
the local magnetization for a spin between two Ji couplings. Accordingly, "asymmetric 
position" denotes a spin position between a Ji and a J2 coupling. 



We have also calculated the local magnetizations for the two models with 
Ji 7^ J2. Local moments can, e.g., be probed with NMR [57, 58] or XMCD 
[59]. The local magnetizations for the Ji = J2 model would simply be pro- 
portional to the total magnetization curve. The results are shown in Fig. 7. 
The calculation demonstrates that the local magnetization of an interacting 
spin system can deviate substantially from the average and even point into 
the opposite direction. 
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